Nonlinear  Optimal  Tracking  Control  of  a 
Piezoelectric  Nanopositioning  Stage 

William  S.  Oates  1  and  Ralph  C.  Smith  2 

Center  for  Research  in  Scientific  Computation 
Department  of  Mathematics 
North  Carolina  State  University 
Raleigh,  NC  27695 

Abstract 

High  performance  nanopositioning  stages,  used  in  a  variety  of  applications  such  as  atomic  force  microscopy 
and  three-dimensional  nanometer-scale  lithography,  require  stringent  position  control  over  relatively  large  dis¬ 
placements  and  a  broad  frequency  range.  Piezoelectric  materials,  which  are  typically  employed  in  nanoposi- 
tioining  stages,  provide  excellent  position  control  when  driven  at  relatively  low  frequency  and  low  field  levels. 
However,  in  applications  where  the  stage  operates  over  a  relatively  large  region  (microns  to  millimeters)  and 
broad  frequency  range  (Hz  -  kHz),  piezoelectric  materials  often  exhibit  nonlinear  and  rate-dependent  hystereis 
which  requires  control  designs  that  can  effectively  accommodate  such  behavior.  In  this  paper,  a  nonlinear, 
thermal-relaxation,  piezoelectric  constitutive  law  is  incorporated  into  an  open  loop  optimal  tracking  control 
design  to  accurately  track  a  desired  reference  signal  when  nonlinearities,  thermal  relaxation  and  hyteresis  are 
present.  A  comparison  between  linear  optimal  control  and  the  nonlinear  optimal  control  design  is  given  to 
illustrate  performance  enhancements  when  the  constitutive  behavior  is  included  in  the  control  design. 
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1.  Introduction 

Highly  accurate  nanopositioning  control  becomes  increasingly  important  in  applications  such  image  steering 
devices  [24],  piezoresponse  force  microscopy  [5],  atomic  manipulation  [22],  and  protein  delivery  [23].  Often, 
piezoelectric  materials  are  chosen  as  actuators  in  nanopositioning  stages  as  well  as  embedded  actuators  within 
MEMS  devices  [24] .  These  materials  provide  small  displacements  and  large  forces  over  a  relatively  large  bandwith 
(Hz  -  MHz);  however,  nonlinearities  and  hysteresis  often  occur  at  moderate  to  large  electric  field  levels.  This 
behavior  requires  the  development  of  a  control  design  that  ensures  adequate  control  input  is  available  to  meet  the 
performance  objectives  while  accommodating  constitutive  nonlinearities  and  hysteresis.  For  example,  the  real¬ 
time  monitoring  of  protein  unfolding  using  an  atomic  force  microscope  requires  highly  accurate  position  control 
at  high  scan  rates  which  are  sufficient  to  introduce  constitutive  nonlinearities  and  hysteresis  in  the  piezoelectric 
actuator  and  subsequently  in  the  positioning  stage.  In  addition,  low  scan  rates  or  static  displacement  control 
can  lead  to  drift  from  relaxation  mechanisms  thus  reducing  accuracy  in  characterizing  a  specimen. 

In  applications  where  high  performance  and  accuracy  are  not  critical,  constitutive  nonlinearities  and  hys¬ 
teresis  can  be  compensated  using  classical  Proportional-Integral  (PI)  or  Proportional-Integral-Derivative  (PID) 
control,  but  this  can  potentially  lead  to  bandwidth  limitation  and  inefficiencies.  Alternatively,  it  is  illustrated  in 
[7,  8]  that  the  use  of  charge-  or  current-controlled  amplifiers  can  essentially  eliminate  hysteresis.  Unfortunately, 
this  mode  of  operation  can  be  prohibitively  expensive  relative  to  more  commonly  employed  voltage-controlled 
amplifiers.  Furthermore,  current  control  is  ineffective  if  maintaining  DC  offsets  as  is  the  case  when  the  x-stage 
of  an  AFM  is  in  a  fixed  position  while  a  sweep  is  performed  with  the  y-stage.  The  present  analysis  focuses  on 
control  development  for  these  applications  where  nanoscale  position  accuracy  is  paramount  and  charge  control 
is  infeasible. 
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Considerable  attention  has  been  focused  on  appropriate  control  strategies  to  compensate  for  ferroelectric 
hysteresis  and  nonlinearities  in  piezoelectric  actuators  used  in  nanopositioning  stages  [3,  12,  13,  14,  19,  25].  PI 
and  PID  control  designs  have  been  shown  to  be  limited  in  bandwidth  and  are  only  effective  in  small  scanning 
regimes  where  hysteresis  is  minimal.  At  moderate  to  large  field  levels  where  the  nonlinearties  and  hysteresis 
increase,  instabilities  may  occur  as  the  control  gains  are  necessarily  increased  to  achieve  the  required  accuracy. 
The  application  of  linear  robust  control  techniques  such  as  TLoo  control  have  provided  improvements  in  bandwidth 
and  robustness  to  nonlinearities  and  hysteresis  [12,  13,  14].  However,  when  linear  robust  control  methods  are 
implemented,  the  control  input  focuses  on  compensating  for  the  ferroelectric  constitutive  behavior  which  reduces 
the  ability  to  compensate  for  external  disturbance  loads  and  system  dynamics  of  the  device.  This  has  led  to 
research  focused  on  the  design  and  implemention  of  nonlinear,  model-based  control  designs  to  accommodate 
constitutive  nonlinearities  and  hysteresis  [3,  11,  19]. 

The  development  and  implementation  of  the  nonlinear  control  design  is  presented  as  follows.  In  Section  2 
a  ferroelectric  homogenized  energy  model  is  summarized  and  employed  in  a  set  of  partial  differential  equations 
(PDE)  and  subsequent  finite  element  model  of  the  piezoelectric  nanopositioning  stage.  In  Section  3,  the  nonlinear 
control  design  is  developed.  First,  open  loop  linear  tracking  control  design  is  developed  and  compared  to  the 
nonlinear  optimal  control  design. 

2.  Nanopositioning  Stage  and  Hysteresis  Model 

The  control  design  is  focused  on  nanopositioning  stages  driven  by  piezoelectric  stack  actuators  for  applica¬ 
tions  in  atomic  force  microscopy  (AFM).  The  typical  AFM  system  is  illustrated  in  Figure  1.  The  side  view 
illustrates  the  nanopositioning  stage  and  piezoelectric  actuator  that  controls  specimen  height  relative  to  the 
cantilever.  In  this  mode  of  operation,  a  photodiode  measures  changes  in  height  of  the  cantilever  position  by 
reflecting  a  laser  beam  off  the  cantilever  and  a  feedback  law  is  used  to  reposition  the  specimen  to  maintain 
constant  surfaces.  The  surface  is  scanned  in  this  manner  which  yields  the  surface  structure  of  the  specimen. 
In-plane  position  control  is  achieved  using  two  additional  stack  actuators  as  illustrated  in  the  top  view  of  the 
positioning  stage  in  Figure  1  (ii) .  In-plane  actuator  coupling  is  typically  negligible  in  this  design  which  allows 
development  of  the  control  design  for  a  single  actuator.  Piezoelectric  tube  actuators  can  also  be  employed  to 
control  in-plane  and  out-of-plane  displacement.  Tube  designs  have  provided  improved  linear  behavior  but  the 
non-negligible  coupling  between  in-plane  displacements  further  complicates  the  control  design  [19].  The  control 
design  presented  here  is  suitable  for  accommodating  coupling  present  in  a  piezoelectric-tube  actuator,  but  the 
finite  element  implementation  is  more  complex;  e.g.,  see  [17].  The  focus  here  is  on  control  development,  therefore 
the  decoupled  stack  actuator  nanopositioner  design  is  considered. 

The  model  development  employed  in  the  present  analysis  focuses  on  a  homogenized  energy  framework  to 
formulate  a  highly  accurate  displacement  control  that  can  be  efficiently  implemented  in  real-time  applications. 
The  constitutive  model  is  based  on  a  previously  developed  homogenized  energy-based  model  [16,  18,  20,  21].  The 
model  incorporates  mesocopic  material  behavior  at  the  domain  or  grain  level  in  a  stochastic  homogenization 
framework  to  predict  macroscopic  material  behavior.  A  distribution  of  interaction  fields  and  coercive  fields 
are  implemented  to  model  polarization  switching  processes  that  typically  occur  in  the  presence  of  material 
inhomogeneities  and  residual  fields.  Boltzmann  relations  are  included  to  model  thermal  relaxation  behavior  when 
thermal  energy  affects  polarization  switching.  Macroscopic  material  behavior  is  determined  by  homogenizing 
the  local  polarization  variants  according  to  the  distribution  of  interaction  and  coercive  fields. 

2.1  Homogenized  Energy  Model 

The  equations  governing  the  homogenized  energy  model  are  summarized  here.  A  detailed  review  of  the  mod¬ 
eling  framework  is  given  in  [16,  18,  20].  The  constitutive  law  is  focused  on  uniaxial  loading  of  rod-type  actuators 
since  the  decoupled  piezoelectric  stack  actuators  are  assumed  in  the  control  design  of  the  nanopositioning  stage 
for  in-plane  position  control.  The  Gibbs  energy  at  the  mesoscopic  length  scale  is 

G(P,  T)  =  ^(P,  T)  -  EP  (1) 

where  ^(P,  T)  is  the  Helmholtz  energy  detailed  in  [16],  T  is  temperature,  E  is  the  electric  field,  and  P  is 
the  polarization.  In  the  one-dimensional  case  considered  here,  the  Helmholtz  energy  function  is  a  double-well 
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Figure  1:  Schematic  of  an  atomic  force  microscope  configuration  used  in  constructing  the  nonlinear  control  design, 
(i)  Side  view  of  the  AFM  set-up  with  out-of-plane  displacement  control  piezoelectric  actuator,  (ii)  Top  view  of  the 
positioning  stage  illustrating  the  configuration  of  the  in-plane  displacement  control  piezoelectric  actuators. 


potential  below  the  Curie  point  Tc  which  gives  rise  to  stable  spontaneous  polarization  with  equal  magnitude  in 
the  positive  and  negative  directions. 

Thermal  relaxation  in  the  piezoelectric  material  is  often  present  and  must  be  addressed  in  the  constitutive 
model  to  predict  creep  that  may  occur  in  the  nanopositioning  stage.  This  can  be  accomplished  by  employing 
the  Boltzmann  relation 


H(G)  =  Ce~GV/kT  (2) 

which  quantifies  the  probability  /r  of  achieving  an  energy  level  G.  The  relative  thermal  energy  ( kT/V )  is 
defined  over  a  representative  volume  element  V  at  the  mesoscopic  length  scale,  k  is  Boltzmann’s  constant,  and 
the  constant  C  is  specified  to  ensure  integration  to  unity.  The  inclusion  of  thermal  energy  in  the  energy  function 
incorporates  the  effect  of  switching  prior  to  the  minima  in  G  disappearing  as  the  temperature  increases.  This 
reduces  the  sharp  transition  of  ferroelectric  switching  near  the  coercive  field  as  the  thermal  energy  increases. 
The  Boltzmann  relation  gives  rise  to  the  expected  values 

r  oo  r  —  Pl 

(P+)  =  /  P^{G)dP  ,  (P_>  =  /  Pn{G)dP  (3) 

J  Pi  J  —  oo 

of  the  polarization  associated  with  positively  and  negatively  oriented  variants.  Here  ±Pj  are  the  positive  and 
negative  inflection  points  in  the  Helmholtz  energy  definition.  Kinetic  equations  are  developed  that  define  the 
volume  fraction  of  positive  x+  and  negative  polarization  variants  x_.  The  kinetic  equations  rely  on  transition 
likelihoods  that  define  the  probability  that  polarization  variants  switch  according  to  the  energy  relations;  see 
[16,  18,  21]  for  details. 

The  resulting  local  average  polarization  is  quantified  by  the  relation 

P  =  x+(P+)  +  x_(P_).  (4) 

The  macroscopic  polarization  is  computed  from  the  distribution  of  local  variants  from  the  relation 

/OO  poo 

/  v(Eo,EI)\P(E  +  EI-,Ec,0](t)dEIdED  (5) 

-oo  J 0 
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where  v{Ec ,  Ei)  denotes  the  distribution  of  coercive  and  interaction  fields  and  £  represents  the  initial  distribution 
of  the  local  variants.  The  local  polarization  P  (E  +  Ej\  Ec,£)  is  determined  from  (4)  where  the  dependence  on 
interaction  fields  and  coercivecl  fields  is  embedded  in  the  kinetic  equations  that  describe  x+  and  x _.  The  density 
i/(Ec,Ei)  is  chosen  as 

v(Ec,  Ej)  =  Cie-h(^/Sc)/2c]2e-£;J2/2b2  (6) 

where  Ec  is  the  average  coercive  field,  c  quantifies  the  coercive  field  variability,  b  is  the  variance  of  the  interaction 
field,  and  c\  is  a  scaling  parameter.  Techniques  to  identify  general  densities  can  be  found  in  [18].  Numerical 
techniques  to  approximate  (5)  and  comparison  to  experimental  results  can  be  found  in  [16]. 

Whereas  the  macroscopic  polarization  is  quantified  by  (5),  the  forces  generated  by  the  stack  actuator  must 
be  quantified  for  implementation  within  the  control  design.  This  is  provided  by  the  constitutive  law 

a  =  Yps  +  cDe  -  hx{P{E)  -  Pr )  -  h2{P{E)  -  Pr )2  (7) 

representing  uniaxial  stress  in  the  piezoelectric  stack  actuator  where  the  effective  properties  of  the  actuator 
include  Yp  as  the  elastic  modulus  at  constant  polarization,  cp  as  the  Kelvin- Voigt  damping  parameter,  e  as  the 
linear  strain  component  in  the  direction  of  loading,  hi  as  the  piezoelectric  coefficient  and  h2  as  the  electrostrictive 
coefficient.  It  is  assumed  that  stress  fields  are  limited  to  the  linear  elastic  regime  where  ferroelastic  switching 
(i.e. ,  stress  induced  domain  switching)  is  negligible.  The  polarization  P(E)  is  computed  using  (5)  where  Pr 
is  the  initial  macroscopic  remanent  state  of  the  material.  In  the  simulations  presented,  the  actuator  model  is 
intially  poled  with  a  field  equal  to  2 Ec  which  results  in  Pr  =  0.206  C/m2. 

Simulations  of  minor  loop  hysteresis  are  illustrated  in  Figure  2  for  the  case  of  zero  applied  stress.  Model 
predictions  illustrate  drift  in  the  polarization  and  strain  behavior  for  an  input  sinusoidal  waveform  at  1  Hz. 
The  material  parameters  associated  with  the  homogenized  energy  model  are  given  in  Table  1.  The  size  of  the 
representative  volume  element  V  was  selected  to  give  typical  relaxation  behavior  for  an  isothermal  process  and 
internal  damping  behavior  at  room  temperature. 

The  stress  computed  using  (7)  includes  linear  stress-strain  behavior  as  well  as  nonlinear  and  hysteretic 
dependence  on  the  electric  field  through  the  P(E)  relation.  It  does  not  include  spatial  dependence.  This  is 
incorporated  in  the  structural  model  in  the  following  section. 

2.2  Structural  Model 

The  constitutive  relations  (5)  and  (7)  are  used  to  develop  a  system  model  that  quantifies  forces  and  dis¬ 
placements  in  the  presence  of  applied  fields  and  stress  that  occur  during  operation  of  the  nanopositioning  stage. 
The  partial  differential  equation  (PDE)  model  is  first  given  and  then  formulated  as  set  of  ordinary  differential 
equations  (ODEs)  through  a  finite  element  discretization  in  space.  The  structural  coupling  of  the  nanoposition¬ 
ing  stage  is  modeled  as  a  damped  oscillator  to  account  for  boundary  conditions  at  the  end  of  the  piezoelectric 
actuator.  The  geometry  used  in  constructing  the  structural  model  is  illustrated  in  Figure  3. 

The  equation  of  motion  for  the  structural  model  is  given  by  the  relation  [4,  16] 


Table  1:  Parameters  employed  in  the  homogenized  energy  model. 

Cl  =  0.54  C/MV2 
c  =  0.4 

b  =  0.87  MV/m. 

Xe  =  0.008  C/mV 
T  =  300  I< 

V  =  1.25  x  105  nm3 
t  =  100  ns 
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(a)  (b) 

Figure  2:  Nonlinear  and  hysteretic  homogenized  energy  model  results  for  minor  loop  hysteresis  typically  observed  in 
ferroelectric  materials  under  zero  stress,  (a)  Electric  field  versus  polarization,  (b)  Electric  field  versus  longitudinal 
microstrain. 


d2w  _  dNtot 
dt2  ~  dx 


(8) 


where  the  mass  density  of  the  actuator  is  denoted  by  p,  the  cross-section  area  is  A  and  the  displacement  is 
denoted  by  w.  The  total  force  Ntot  acting  on  the  actuator  is 


,  dw 


d2 


w 


Ntot(t,  x)  =  Y^A—  +  cdA—  +  Fp(E)  +  Fd 


dx 


dxdt 


(9) 


where  the  the  first  term  on  the  right  hand  side  of  the  equation  represents  the  elastic  restoring  force  and  a  linear 
damping  coefficient  is  incorporated  in  the  second  term.  The  linear  elastic  strain  component  in  the  direction  of 
loading  is  defined  by  e  =  The  term  Fd  incorporates  external  disturbance  loads  and  the  coupling  force  Fp 
represents  forces  generated  by  an  applied  electric  field  where 


+w 


x=0  d  x=L 


kLw 

- 


cL  dw/dt 


/ 

/ 

/ 

/ 

/ 

/ 

/ 

V 


Figure  3:  Piezoelectric  stack  actuator  with  damped  oscillator  used  to  quantify  loads  during  scanning  operations  of  an 
AFM  nanopositioning  stage.  Disturbance  forces  along  the  actuator  are  given  by  Fd  and  the  control  input  is  u(t). 
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Table  2:  Model  parameters  for  the  piezoelectric  stack  actuator  and  damped  oscillator. 


Yp  =  60  x  109  N/m2 
h i  =  1.0  x  106  N/C 
h2  =  1.0  x  103  Nm2/C' 
PR  =  0.3  C/m 2 


p=  7.5  x  103  kg/m3 
A  =  5.07  x  10"4  m2 
cd  =  3.7  x  106  Ns/m 
L  =  0.1  to 


kL  =  3.04  x  106  N/m 
cl  =  3.04  x  102  Ns/m 
mL  =40  g 


FP(E )  =  A[hi(P(E)  -  Pr)  +  h2(P(E )  -  Pr )2]  (10) 

and  the  hysteretic  and  nonlinear  E  —  P  relation  is  specified  by  (5). 

As  illustrated  in  Figure  3,  the  boundary  conditions  are  defined  by  a  zero  displacement  at  x  =  0  and  the 
balance  of  forces  between  the  piezoelectric  actuator  and  the  damped  oscillator  at  x  =  L  yields 

dw  d2w 

Ntot(t,L)  =  - kLw(t,L )  -  cL  —  (t,  L)  -  mL-^p(t,L).  (11) 

The  initial  conditions  are  w(0,x)  =  0  and  =  0.  Model  parameters  associated  with  the  stack  actuator 

and  damped  oscillator  used  in  the  control  design  are  given  in  Table  2. 

The  strong  form  of  the  PDE  model  given  by  (8)  can  be  written  in  the  weak  form  for  finite  element  imple¬ 
mentation;  see  [16]  for  details.  The  weak  form  of  the  model  given  is  used  to  obtain  a  matrix  ODE  system  given 

by 


Mw  +  C^w  +  Kw  =  Fp(E)  b  +  fj  (12) 

where  w (<)  =  [w\(t), . . .  ,WN(t)]  are  the  nodal  solutions  obtained  from  the  weak  formulation  and  finite  element 
discretization.  M  G  RNxN ,  Cr  €  RNxN  and  K  G  M.NxN  denote  the  mass,  damping  and  stiffness  matrices.  The 
vectors  b  G  and  fj  G  include  the  integrated  basis  functions  related  to  the  control  input  and  disturbance 
loads,  respectively.  Components  contained  within  the  system  matrices  and  vectors  can  be  found  in  [11,  16]. 
Formulation  of  (12)  as  a  first  order  system  yields 

x(f)  =  Ax(f)  +  [B(w)](t)  +  G(t) 

x(0)  =  x0  (13) 

y(t)  =  Cx(t) 


where  x(f)  =  [w,  w]T  should  not  to  be  confused  with  the  coordinate  x.  The  matrix  A  incorporates  the  mass, 
damping  and  stiffness  matrices  given  in  (12)  and  [B(u)](t)  includes  the  nonlinear  input  where  u(t)  is  defined  as 
the  electric  field.  The  initial  conditions  are  defined  by  x0.  The  output  of  the  system  y(t)  is  a  function  of  the 
system  states  according  to  the  matrix  C.  In  the  nanopositioning  stage,  it  is  assumed  that  only  the  displacement 
at  x  =  L  is  observable  which  results  in  C  =  [  1  0  ...  0  ]  with  dimension  1  x  2 N .  External  disturbances 

are  incorporated  in  G. 

A  temporal  discretization  of  the  system  given  by  (13)  is  used  to  numerically  analyze  the  dynamic  perfor¬ 
mance.  The  trapezoid  rule  is  adopted  since  it  is  moderately  accurate,  A-stable,  and  requires  minimal  computer 
storage  capacity.  The  trapezoidal  discretization  yields  the  relation  for  each  iteration 


xfc+1  =  Wxfc  +  V  [B(ufc)]  +  [G (4)  +  G(ifc+1)] 
x(0)  =  x0 


(14) 
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where  a  temporal  step  size  At  is  employed  giving  a  discretization  in  time  defined  by  tk  =  kAt.  The  values  x*, 
approximate  x(tfc).  The  matrices 


W  = 


V  =  At 


-l 


(15) 


are  created  once  for  numerical  implementation  yielding  approximate  solutions  with  0(h2,  (At)2)  accuracy.  In 
(14),  only  values  for  the  control  input  at  the  present  time  tk  are  included  in  the  temporal  discretization  to  model 
dynamics  of  real-time  feedback  developed  in  the  following  section.  This  approach  is  employed  since  the  control 
input  can  only  be  dependent  on  the  current  and  previous  states. 

3.  Control  Design 

We  compare  the  performance  between  linear  and  nonlinear  optimal  control  designs  when  ferroelectric  non- 
linearities,  hysteresis  and  relaxation  behavior  are  present  in  the  piezoelectric  actuator.  The  ferroelectric  consti¬ 
tutive  behavior  is  shown  to  introduce  significant  degradation  in  position  accuracy  when  linear  optimal  control 
methods  are  employed.  Improvements  in  high  accuracy  tracking  are  then  obtained  by  implementing  a  nonlinear 
optimal  control  design.  Actuator  displacements  on  the  order  of  tens  of  microns  are  chosen  for  the  analysis 
where  it  has  been  demonstrated  that  nonlinearities  and  hysteresis  are  significant.  The  reference  tracking  profile 
used  in  the  simulations  is  based  on  a  typical  scan  and  hold  procedure  used  in  nano-  and  micron-scale  materal 
characterization.  This  procedure  requires  that  the  specimen  translate  in  the  x-direction  while  holding  in  the 
y-direction.  At  the  end  of  each  scan,  the  position  is  held  in  the  x-direction  and  increased  or  decreased  a  small 
amount  in  the  y-direction  and  the  procedure  is  repeated  in  the  negative  ^-direction.  A  low  frequency  scan  rate 
is  chose  for  the  simulations  (0.2  Hz)  where  relaxation  and  creep  are  significant. 


3.1  Optimal  Tracking  Control 

Development  of  the  nonlinear  optimal  tracking  control  design  follows  a  previous  approach  focused  on  negli¬ 
gible  relaxation  hysteresis  of  magnetostrictive  actuators  for  vibration  attenuation  of  beam  and  plate  structures 
and  tracking  control  of  rod  structures  [10,  11,  15].  We  summarize  here  key  equations  associated  with  opti¬ 
mal  tracking  control  and  its  application  in  compensating  ferroelectric  nonlinearities,  hysteresis  and  relaxation 
behavior. 

Optimal  tracking  control  utilizes  a  cost  functional  to  determine  the  optimal  control  input.  The  cost  functional 


J  =  ^(Cx(^/)  -  r(tf))T P{Cx{tf)  -  r(tf))  + 


lS  r 


H  -  X T(t)k(t) 


dt 


(16) 


penalizes  the  control  input  and  the  error  between  the  piezoelectric  actuator  displacement  and  the  prescibed 
displacement  where  P  penalizes  large  terminal  values  on  the  tracking  error,  H  is  the  Hamiltonian,  and  A (£)  is 
a  set  of  Lagrange  multipliers.  The  Hamiltonian  is 


H  =  i  [(Cx(t)  -  r(t))TQ(Cx(t)  -  r(t))  +  uT {t)Ru(t)\ 
+XT  [Ax(t)+  |B(u)](i)  +  G(i)] 


(17) 


where  penalties  on  the  tracking  error  and  the  control  input  are  given  by  the  variables  Q  and  i?,  respectively. 

The  minimum  of  the  cost  functional  in  (16)  is  determined  under  the  constraint  of  the  differential  equation 
given  by  (13).  By  employing  Lagrange  multipliers  an  unconstrained  minimization  problem  is  constructed  where 
the  stationary  condition  for  the  Hamiltonian  yields  the  adjoint  relation  [2,  6] 
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A (t)  =  -AtA (t)  -  CTQCx(t)  +  CTQr(t) 


(18) 


and  optimal  control  input 


The  resulting  optimality  system  is 


'x(f)  ' 

Ax(t)  +  [B(u)](t)  +  G  (t) 

.  Ht)  . 

-AtA (t)  -  CTQCx(t)  +  CTQr(t)  _ 

x(t0)  =  x0 

\(tf)  =  CTP(Cx(tf)-r(tf)). 


(19) 


(20) 


The  force  determined  from  (10)  is  included  in  the  input  operator  [B(u)](f)  which  directly  includes  the  nonlinear 
and  hysteretic  E  —  P  relation  as  well  as  relaxation  behavior  within  the  control  formulation.  This  system  of 
equations  results  in  a  two-point  boundary  value  problem  which  presents  challenges  in  obtaining  a  solution  for 
large  systems.  Furthermore,  the  nonlinear  nature  of  the  input  operator  precludes  an  efficient  Riccati  formulation. 


3.2.1  Linear  Optimal  Control 

If  the  input  operator  is  assumed  to  be  linear  and  not  a  function  of  the  polarization  state  (i.e.  valid  for  small 
to  moderate  field  inputs) ,  the  control  design  can  be  simplified  by  determining  the  suboptimal  control  input  from 
the  algebraic  Riccati  equation 


ATn  +  lea  nBR-1BTn  +  ctqc  =  o 


see  [2,  6]  for  details. 

The  linear  optimal  control  input  is  then  defined  by 


u*{t)  =  —R  1BT[IIx(i)  —  v(t)\ 


where  the  variable  h>(t)  £  R2Ar  is  the  solution  to  the  auxiliary  differential  equation 

i >(t)  =  -  [A  -  BiT^n] T  i !/(t)  -  C TQr(t)  +  IIG(i) 
u{tf)  =  CTP r(tf) . 


(21) 


(22) 


(23) 


3.2.2  Nonlinear  Open  Loop  Optimal  Tracking  Control 

When  nonlinearities  in  the  input  operator  cannot  be  neglected,  (20)  has  the  general  form 

z(f)  =  F(t,  z) 

E0z(t0)  =  [x0,0]T  (24) 

Efz(t/)  =  [0,0]T 


where  z  =  [x(t),A(t)]T  and 


(25) 


A  x(t)  +  [B(u)](f)  +  G  (t) 

-ATX (t)  -  CTQCx(t)  +  C TQr{t) 


Eq  — 


I  0 
0  0 


Ef 


0  0 
-CTPC  I 


Here  I  denotes  an  identity  matrix  with  dimension  corresponding  to  the  number  of  basis  functions  employed  in 
the  spatial  approximation  of  the  state  variables.  Also  note  that  the  prescribed  displacement  has  been  restricted 
to  r(tf)  =  0. 

The  system  given  by  (24)  is  approximated  by  discretizing  the  time  interval  [to,tf]  with  a  uniform  step  size 
At  at  the  points  , tjv  =  tf.  The  approximate  values  of  the  state  solutions  and  the  adjoint  at  each 

time  step  are  denoted  by  Zo,  •  •  •  ,  zn-  Whereas  several  techniques  are  available  to  approximate  the  solution  to 
(24),  such  as  finite  differences  and  nonlinear  multiple  shooting  [1],  a  central  difference  of  the  temporal  derivative 
is  utilized  here  which  provides  a  formulation  that  can  be  cast  in  a  analytic  LU  decomposition  to  solve  the 
two-point  boundary  value  problem,  details  describing  the  numerical  approach  are  described  elsewhere  [11].  This 
approach  results  in  expressing  (24)  as  the  problem  of  finding  Zh  =  [zo,  •  •  •  ,  z n]  which  solves 

^(zh)  =  0.  (26) 

Equation  (26)  includes  the  optimality  system  at  each  time  step  and  the  boundary  conditions  given  in  (25). 
Details  are  given  in  [15]. 

A  quasi-Newton  iteration  of  the  form 


where  solves 


(27) 

(28) 

is  then  used  to  approximate  the  solution  to  the  nonlinear  system  given  by  (26). 


3.2.3  Simulation  Results 

Improvements  in  tracking  control  are  demonstrated  here  by  comparing  the  linear  and  nonlinear  optimal 
control  designs.  The  values  used  to  penalized  the  displacement  error  and  control  input  were  Q  =  lx  1016 
and  R  =  lx  10-8.  The  penalty  on  the  final  state  was  P  =  1.4  x  1010.  In  Figure  4,  linear  optimal  control 
illustrates  poor  tracking  performance  using  a  scan  rate  of  0.2  Hz.  The  thermal  relaxation  effect  is  significant  in 
regions  of  static  displacement  as  seen  in  Figure  4(a).  In  Figure  5,  the  nonlinear  control  design  is  applied  to  the 
piezoelectric  actuator  again  using  a  scan  rate  of  0.2  Hz.  Significant  improvements  are  achieved  in  tracking  the 
reference  signal.  In  Figure  5(a)  the  commanded  displacement  and  error  between  the  commanded  displacement 
and  reference  displacement  are  given.  Whereas  relaxation  behavior  is  present  at  this  operating  frequency  as 
illustrated  in  Figure  5(b),  the  nonlinear  control  design  significantly  reduces  error  relative  to  the  linear  control 
design.  This  is  clearly  demonstrated  in  regions  of  static  displacement  where  the  control  input  increases  or 
decreases  to  ensure  the  polarization  is  constant  which  gives  rise  to  static  displacement. 

4.  Concluding  Remarks 

An  open  loop  nonlinear  optimal  control  design  has  been  applied  to  a  nanopositioning  stage  for  quasi-static 
scanning  rates  for  applications  where  nanoscale  resolution  is  important.  A  homogenized  energy  model  was 
implemented  in  the  control  design  to  account  for  local  residual  interaction  fields  and  variations  in  the  coercive 
field.  Relaxation  behavior  was  incorporated  into  the  constitutive  model  to  account  for  creep  during  quasi-static 
actuation.  The  nonlinear  control  design  significantly  improved  reduction  in  tracking  errors  relative  to  linear 
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Figure  4:  Tracking  performance  using  linear  optimal  control  to  compensate  for  relaxation  behavior  at  0.2  Hz.  (a) 
Commanded  displacement  and  the  desired  reference  signal,  (b)  The  associated  control  input  and  polarization  relaxation 
behavior. 


optimal  control.  However,  open  loop  control  can  lead  to  deleterious  effects  when  operating  uncertainties  are 
present.  Improved  robustness  can  be  achieved  by  introducing  perturbation  feedback  around  the  optimal  open 
loop  control  input  but  this  is  beyond  the  scope  of  the  paper.  Details  decribing  performance  robustness  for 
tracking  control  using  perturbation  feedback  can  be  found  in  [9,  11]. 


(a) 


Figure  5:  High  performance  tracking  using  nonlinear  optimal  control  to  compensate  for  relaxation  behavior  at  0.2  Hz. 
(a)  Commanded  displacement  and  error  between  the  actuator  response  and  the  reference  signal,  (b)  The  associated 
control  input  and  polarization  relaxation  behavior. 
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